Rar

Among Us Drip. I thought it was funny, idk.

Task 1

Working Directory

getwd()
## [1] "/home/joshua/Desktop/MainFolder/OuClasses/Spring 2023/Applied Statistical Methods/FALL224753wise0046/LABS/LAB5"

Task 2

mybin function

mybin=function(iter=100,n=10, p=0.5){ 
  # make a matrix to hold the samples initially filled with NA's
  sam.mat=matrix(NA,nrow=n,ncol=iter, byrow=TRUE)
  
  #Make a vector to hold the number of successes in each trial
  succ=c()
  
  for(i in 1:iter){
    #Fill each column with a new sample
    sam.mat[,i]=sample(c(1,0),n,replace=TRUE, prob=c(p,1-p))
    
    #Calculate a statistic from the sample (this case it is the sum)
    succ[i]=sum(sam.mat[,i])
  }
  #Make a table of successes
  succ.tab=table(factor(succ,levels=0:n))
  
  #Make a barplot of the proportions
  iter.lab = paste0("iter = ", iter)
  n.lab = paste0("n = ",n)
  p.lab = paste0("p = ",p)
  lab = paste0(iter.lab, n.lab, p.lab, sep=", ")
  barplot(succ.tab/(iter), col=rainbow(n+1), main="Binomial simulation", xlab="Number of successes")
  succ.tab/iter
}

100 iterations Plot

bin100 = mybin(iter = 100, n = 10, p = 0.7)

200 iterations Plot

bin200 = mybin(iter = 200, n = 10, p = 0.7)

500 iterations Plot

bin500 = mybin(iter = 500, n = 10, p = 0.7)

1000 iterations Plot

bin1000 = mybin(iter = 1000, n = 10, p = 0.7)

10000 iterations Plot

bin10000 = mybin(iter = 10000, n = 10, p = 0.7)

10000 iteration table

names(bin10000) = 0:10
bin10000
##      0      1      2      3      4      5      6      7      8      9     10 
## 0.0000 0.0000 0.0014 0.0104 0.0379 0.1103 0.1999 0.2654 0.2265 0.1187 0.0295

Verify Values

dbin.tab = round(dbinom(0:10, size = 10, prob = 0.7), 4)
names(dbin.tab) = 0:10
dbin.tab
##      0      1      2      3      4      5      6      7      8      9     10 
## 0.0000 0.0001 0.0014 0.0090 0.0368 0.1029 0.2001 0.2668 0.2335 0.1211 0.0282

Task 3

Sample without Replacement

marbles = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0)
sample(marbles, size = 5, replace = FALSE, prob = NULL)
## [1] 1 0 1 1 0

Sample with Replacement

sample(marbles, size = 5, replace = TRUE, prob = NULL)
## [1] 1 0 0 1 1

myhyper function

myhyper=function(iter=100,N=20,r=12,n=5){
  # make a matrix to hold the samples initially filled with NA's
  sam.mat=matrix(NA,nr=n,nc=iter, byrow=TRUE)
  #Make a vector to hold the number of successes over the trials
  succ=c()
  for( i in 1:iter){
  #Fill each column with a new sample
  sam.mat[,i]=sample(rep(c(1,0),c(r,N-r)),n,replace=FALSE)
  #Calculate a statistic from the sample (this case it is the sum)
  succ[i]=sum(sam.mat[,i])
  }
  #Make a table of successes
  succ.tab=table(factor(succ,levels=0:n))
  #Make a barplot of the proportions
  barplot(succ.tab/(iter), col=rainbow(n+1), main="HYPERGEOMETRIC simulation", xlab="Number of successes")
  succ.tab/iter
}

100 iterations Plot

hyper100 = myhyper(iter=100, n=5, N=20, r=12)

200 iterations Plot

hyper200 = myhyper(iter=200,n=5, N=20,r=12)

500 iterations Plot

hyper500 = myhyper(iter=500,n=5, N=20,r=12)

1000 iterations Plot

hyper1000 = myhyper(iter=1000,n=5, N=20,r=12)

10000 iterations Plot

hyper10000 = myhyper(iter=10000,n=5, N=20,r=12)

10000 iteration table

names(hyper10000) = 0:5
hyper10000
##      0      1      2      3      4      5 
## 0.0041 0.0508 0.2366 0.4018 0.2562 0.0505

Verify Values

dhype = round(dhyper(x=0:5, m=12, n=8, k=5), 4)
names(dhype) = 0:5
dhype
##      0      1      2      3      4      5 
## 0.0036 0.0542 0.2384 0.3973 0.2554 0.0511

Task 4

What does mysample() do?

mysample=function(n, iter=10,time=0.5){
  for( i in 1:iter){
    #make a sample
    s=sample(1:10,n,replace=TRUE)
    # turn the sample into a factor
    sf=factor(s,levels=1:10)
    #make a barplot
    barplot(table(sf)/n,beside=TRUE,col=rainbow(10), 
    main=paste("Example sample()", " iteration ", i, " n= ", n,sep="") ,
    ylim=c(0,0.2)
    )
    
    #release the table
    Sys.sleep(time)
  }
}

The mysample function takes a vector n, number of iterations iter, and amount of time time. Basically it will sample 10 random samples from the vector, turn that sample to a factor, and then make barplot of that factor. Then the function will wait time time seconds and then loop iter-1 more times after the inital iteration.

Run Line

# Original Line
# mysample(n=1000, iter=30, time=1)
# Changed Line
mysample(n=1000, iter=1, time=0)

The line takes 30 seconds to run because it waits 1 second after each iteration, in this case 30 (30*1 = 30 seconds). But it also samples 10 random samples from the 1000 sided population vector.

Task 5

R Calculations

Calculation 1

choose(8,4)
## [1] 70

Calculation 2

1-ppois(4, lambda = 2)
## [1] 0.05265302

More R Calculations

More Calculations 1

# dnbinom(x=10-3,size=3,prob=0.4)
mynbin=function(y,r,p){
choose(y-1,r-1)*p^r*(1-p)^(y-r)
}
mynbin(10,3,0.4)
## [1] 0.06449725

More Calculations 2

pbinom(8, 15, 0.4)
## [1] 0.9049526

Task 6

Load Custom Package

bin = BestRPackage::abin(100, 10, 0.5)

bin
## 
##    0    1    2    3    4    5    6    7    8    9   10 
## 0.00 0.00 0.07 0.20 0.23 0.17 0.14 0.16 0.03 0.00 0.00

Task 7

negBinomialDistribution=function(iter=100,n=10, p=0.5){ 
  # make a matrix to hold the samples initially filled with NA's
  sam.mat=matrix(NA,nrow=n,ncol=iter, byrow=TRUE)
  
  #Make a vector to hold the number of successes in each trial
  succ=c()
  
  for(i in 1:iter){
    #Fill each column with a new sample
    sam.mat[,i]=sample(c(1,0),n,replace=TRUE, prob=c(p,1-p))
    
    #Calculate a statistic from the sample (this case it is the sum)
    succ[i]=sum(sam.mat[,i])
  }
  #Make a table of successes
  fail.tab=table(factor(n-succ,levels=0:n))
  
  #Make a barplot of the proportions
  iter.lab = paste0("iter = ", iter)
  n.lab = paste0("n = ",n)
  p.lab = paste0("p = ",p)
  lab = paste0(iter.lab, n.lab, p.lab, sep=", ")
  barplot(fail.tab/(iter), col=rainbow(n+1), main="Neg Binomial simulation", xlab="Number of successes")
  fail.tab/iter
}
negBinomialDistribution(100, 10, 0.7)

## 
##    0    1    2    3    4    5    6    7    8    9   10 
## 0.01 0.10 0.29 0.24 0.22 0.12 0.02 0.00 0.00 0.00 0.00

I think this is right, if not, at least I tried.